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Abstract 

A numerical algorithm has been developed for solving the equations 
describing chemically reacting supersonic flows. The algorithm employs a two- 
stage Runge-Kutta method for integrating the equations in time and a Chebyshev 
spectral method for integrating the equations in space. The accuracy and 
efficiency of the technique have been assessed by comparison with an existing 
implicit finite-difference procedure for modeling chemically reacting flows. 
The comparison showed that the new procedure yielded equivalent accuracy on 
much coarser grids as compared to the finite-difference procedure with 
resultant significant gains in computational efficiency. 
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Nomenclature 


cross-sectional area, constant in Arrhenius law 
constants in specific heat equations 
concentration of species 
speed of sound 

specific heat at constant pressure 
activation energy 
total internal energy 
flux vector 

expansion coefficients in Chebyshev series 
mass fraction 
source vector 
total enthalpy 

reference enthalpy at standard conditions 

identity matrix 

equilibrium constant 

reverse reaction rate 

forward reaction rate 

molecular weight 

number of nodes 

number of reactions 

number of species 

static pressure 

total pressure 

steady-state residual 

universal gas constant 

total temperature 

static temperature 

Chebyshev polynomial 

time 

time step 

species production rate 
dependent variable vector 
velocity 



X 


spatial variable 
Ax spatial step size 

Y stoichiometric coefficient, ratio of specific heats 

A eigenvalue 

p density 

<f> equivalence ratio 

Subscripts 

c based on chemistry 

f based on fluids 

i>j species indices 

R reactions, reference value 

s species 

sp evaluated spectrally 

Superscript 

- mass weighted value 


iv 



INTRODUCTION 


Research to develop ramjet and supersonic combustion ramjet (scramjet) 
propulsion systems has been underway at the NASA Langley Research Center for a 
number of years. A critical element in the design of scramjets and ramjets is 
the detailed understanding of the complex flow field present in the engine 
over a range of operating conditions. Numerical modeling of various regions 
of the engine flow field has proven to be a valuable tool for gaining insight 
into the nature of these flows. In recent years, computer programs have been 
developed to model the chemically reacting flow fields in ramjet and scramjet 
systems. 1 ^ ^ These programs have employed both explicit and implicit finite- 
difference procedures to solve the equations governing the flow field in the 
engine combustors. The calculations have often required long computer runs to 
reach desired steady-state conditions and have been quite costly, due to 
stiffness introduced in the equations by the finite-rate chemical kinetics 
that is required for accurate modeling. Also, computer resource limitations 
have sometimes reduced the degree of spatial resolution that could be achieved 
in the calculations. These factors have led to the desire for more efficient 
and more highly accurate algorithms for solving chemically reacting combustor 
flows • 

The system of partial differential equations describing chemically 
reacting flows (see eq. (4)) is stiff because of the highly disparate time 
scales that exist among the equations. Certain chemical reactions in an 

overall combustor kinetics system can take place on an extremely short scale 

-12 -3 
of the order of 10 seconds, whereas the fluids dynamics may require 10 to 

10 seconds for a typical case to reach steady-state conditions. There are of 

course several intermediate scales lying between these two extremes. Mathe- 
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matically, stiffness is often defined by examining the eigenvalues of the 
Jacobian of the governing equation system and noting that the ratio of the 
real part of the largest to real part of the smallest eigenvalue is a large 
number. The former physical definition is perhaps the more useful test of 
stiffness; it is felt directly in the numerical integration of stiff systems 
through the required proper choice of the integration time step. We will deal 
with this requirement now and then follow with a discussion concerning inte- 
gration of the spatial part of the problem. 

Stiffness in the system of equations governing chemically reacting flows 
typically arises from the source terms in the equations describing production 
and loss of the chemical species that are present. Large values for these 
source terms produce rapid changes in the dependent variables being sought, 
and result in the very short time scales discussed in the previous 
paragraph. To explore the problem of mixed (short and long) time scales, 
consider the ordinary differential equation (ODE) system 4 

||-W? (1) 


where f = [fpf 2 ]^> f(0) = [2,1] T and 


A = 


1]^ and 


’ -500.5 

499.5 

499.5 

-500.5 

= -1000.0 

and 


equation (1) follows as 
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r -i c -t , n -lOOOt 

f^(t) = 1.5 e + 0.5 e 


, , -t n , -lOOOt 

f 2 (t) = 1.5 e - 0.5 e 


( 2 ) 


Note that the solutions and have a rapidly decaying component 
corresponding to X^ and a much more slowly decaying component corresponding 
to X^. If we were solving this problem numerically, accuracy would require 
that we advance the solution from the initial conditions using very small time 
steps. However, once the solution dominated by X^ decays, we would prefer to 
advance the solution using larger time steps that would still maintain an 
acceptable level of accuracy. Care must be taken in picking a numerical algo- 
rithm that will allow this choice of time step. Otherwise, the numerical 
stability of the solution will continue to be dictated by X^ even though its 
component has decayed, and very small time steps will still be required to 
maintain stability. In response to this difficulty, several authors, 
including Bussing and Murraan, 5 Stalnaker, et al., 8 and Smoot, Hecker, and 
Williams 7 recognized that the stiff source terms in the system of equations 
governing chemically reacting flow should be evaluated implicitly. Therefore, 
for our problem, algorithms should be developed with the source terms written 
implicitly at the new time level in the integration step. Other terms in the 
governing equations, that do not lead to stiffness, can still be evaluated 
explicitly. 5 5 7 

Perhaps the best known algorithms for solving stiff systems of ordinary 
differential equations are those developed by Gear. 8 These schemes employ 
Adams^ methods of variable order with explicit formulas used to solve non- 
stiff equations in the system and implicit formulas used to solve the stiff 
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equations. Hindmarsh generalized the Gear algorithms ^ and developed a variant 
to allow for a variable integration step size , 10 Another class of algorithms 
for effectively solving stiff ordinary differential equations is the 
exponential-fitted schemes, described in the work of Liniger and 
Willoughby 11 and Pratt . 12 These methods fit the solution at two or more inte- 
gration points with an exponential interpolant rather than polynomial 
interpolants used in the Gear schemes. The exponential-fitted schemes more 
naturally follow the exponential behavior of solutions to chemical kinetics 
problems. Additionally, for decaying solutions, this class of algorithms has 
an infinite stability radius for both the implicit and explicit variants . 12 A 
good deal more work has been carried out to develop efficient and accurate 
algorithms for solving stiff systems of ordinary differential equations 
resulting from chemical kinetics problems. We will not continue our survey 
here, but rather we will refer readers to an interesting paper by Bui, 
Oppenheira, and Pratt 13 that further discusses the area of stiff ODE solvers. 

Next, we deal with the computation of spatial derivatives in the 
governing equations. The importance of accurately modeling spatial deriva- 
tives cannot be overemphasized. Chemical reaction does not take place until 
fuel and oxidant are brought together and macroscopically mixed by convective 
transport, and then mixed down to the microscopic (molecular) level by 
diffusive processes. To model these processes, spatial derivatives must be 
accurately computed. Due to computer storage limitations, higher order 
numerical methods were indicated. Higher order finite-difference schemes 
offered one option for computing the spatial derivative. Another option was 
apparent from earlier work of the second and third authors to develop methods 
for highly accurate solutions of the Euler equations. Here, Hussaini, 
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et.al. 14 used a spectral collocation method to compute the required spatial 
derivatives in the governing equations. Using this approach, several problems 
governed by the Euler equations were successfully solved and accurate 
solutions were obtained on relatively coarse grids as compared to finite- 
difference solutions of the same problems, c Spectral methods) are based on the 
representation of the solution to a problem f by a finite series of global 
functions X of the form 

N „ 

f(x) = I a X (x), (3) 
n=0 n n 

A 

where the are the expansion coefficients of the series.^ The X n should be 
a complete orthogonal set. Spatial derivatives of f are then approximated by 
taking derivatives of the corresponding series (3). If properly applied, the 
high order approximation (3) yields a very accurate numerical representation 
for derivatives of f. Spectral methods should therefore satisfy our 
requirements for approximating spatial derivatives in the equation governing a 
chemically reacting flow field. 

This paper discusses the development of a numerical algorithm for solving 
the equations governing a chemically reacting flow. The algorithm employs a 
two-stage partial implicit Runge-Kutta scheme for integrating the equations in 
time and a Chebyshev spectral collocation method for computing spatial 
derivatives in the equations. A computer program has been written to apply 
this algorithm for the solution of a reacting flow problem. The code is 
currently limited to quasi-one-dimensional inviscid flows with hydrogen-air 
reaction, which is appropriate for development and evaluation of the 
algorithm. There appear to be no restrictions prohibiting extension of the 
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algorithra to three-dimensional viscous flows. Chemical reaction is 
represented in the program with a finite-rate chemistry model, and a real gas 
thermodynamic model is employed. 


ANALYSIS 

Governing Equations 

The quasi-one-dimensional Euler equations in conservation law form with 
multiple species undergoing chemical reaction are 16 




(4) 


where 


-v T 

U = {pA, puA, pe Q A, pf ± a} 

F = {puA, pu 2 A + pA, puh Q A, puf^ A} 


S = {0, -p 


dA(x) 

dx 


, 0, -w ± A} T , 


(5) 

C6) 

(7) 


and 


N 


h = / cdT+~+). 

0 i. p 2 A 


(H ?)j £ j 



( 8 ) 

( 9 ) 


If there are N g chemical species, then, i = 1,2,..., (N g - 1) and (N g - 1) 
equations must be solved for the species f^. The final species mass fraction 
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f N can then be found by conservation of mass since 


N 

s 



Chemistry Model 

The chemical reaction of hydrogen and oxygen is modeled in this work with 
the global finite-rate hydrogen-air chemistry model of Rogers and Chinitz. 17 
This model adequately represents the chemical reaction taking place in the 
problems to be considered, and it also produces an extremely large disparity 
in the time scales present in the problems. This phenomenon allows the 

ability of the numerical algorithm to deal with resulting stiffness to be 
demonstrated. 

The Rogers-Chinitz model assumes that the overall reaction of hydrogen 
and oxygen takes place through two reactions, the first resulting in the 
formation of hydroxyl radical, and the second combining the hydroxyl radical 
with hydrogen to form water. The reactions are given by, 


H 2 + °2 


fl 


bl 


20H 


20H + H 


f 2 

9 2H,0 

z k ^ 

d2 


( 10 ) 


(ID 


where the kf's are the forward reaction rates and the are the reverse 

reaction rates. The reverse rates can be found given the forward rate and 
equilibrium constant K for each reaction, as 
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k, = k/K. 
b r 


( 12 ) 


The forward reaction rates are computed from the Arrhenius Law, 


N . -E . /R°T 
k = A.T 1 e X 
i 1 


(13) 


for each reaction i. For the Rogers-Chinitz model, the rates are given by, 17 


k = A T -10 e " 4865/R ° T 

i 1 A 1 1 6 


(14) 


where 


k = A T" 13 e “42500/R°T 

Kf 2 a 2 1 e 


A. = (8.917 4 + 31,433 ~ 28.95)(10 47 ) Cm3 

1 


(15) 


mole-s 


A 2 - (2.0 + iip- 0.833 <t>) ( 10 64 ) ^r £ 


and 


K, 


= 26.164 e 


-8992/T 


K 


2 


= 2.682 x 10 


Xe 694l5/T^ 


Knowing the reaction rates for reactions (10) and (11), the production of the 
four species present in the model can be found from the law of mass action. 
For a general reaction, 


N 

s 


l 

j=l 


Y ij C j 


k fi N s 

^ 1 Y ijV 

k bi J =1 ' 


i = 1,2,. 


■>K 
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the law of mass action states that the rate of change of concentration of 
species j by reaction i is given by, 18 




N 


N 


ij _ 


k*. E c, iJ - k, , n c; 
. £1 J-l 3 bl j-1 3 


(16) 


The rate change in concentration of species j by all N R reactions is then 
found by summing the contributions from each reaction, 


N R 

C i = l (C.).. (17) 

J i=l J 


Finally, the production rate of species j is found from, 


w. = C.M. . 
J J J 


(18) 


Applying the law of mass action to the global model, equations (10) and (11), 
gives , 17 


^0 2 “ " k fl C H 2 C 0 2 + k bl C 0H * 


S 2 0 2(k f2 C 0H C H 2 “ k b2 C H 2 0 ^ ’ 


^H 2 L 0 2 2 C H 2 0’ 


£ oh (2 £ 0 2 + 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 


The source terms for the last i equations in (4) can now be determined, as a 
function of the dependent variables, by application of equation (18). 



- 10 - 


Thermodypamlcs Model 

The specific heat at constant pressure, c p , is nearly a linear function 
of temperature for each species present in the flow field (H 2 , O 2 , OH, H 2 O, 
N 2 ) over the range of temperature being considered. We therefore fit the c p 
versus temperature data 19 for each species i with 


c (T) = a T + b . , (23) 

V t i 1 


where a and b are constants. A mixture specific heat, c , can then be 
defined by weighting over the species i as 


c 

P 


N 

s 



f . . 

1 


(24) 


The total enthalpy of the mixture, made up of the five species, is given by 


H = 


N 

s 





dT + H 


o 

T. 


1 



(25) 


19 

where H° is the reference enthalpy at the reference temperature T R = 0 K. 
Putting equation (24) into (25) and integrating gives 


H 


N 

s 


l 

i=l 



+ b T + H° 

i 



(26) 


Finally, the mixture gas constant, R, is found by weighting the individual gas 
constants over the species i as 


N 

s 


i=l 


R i f i' 


R = 


(27) 
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Equations (24), (26), and (27) can then be used to define all other required 
thermodynamic variables. 


SOLUTION OF THE GOVERNING EQUATIONS 

Chebyshev Spectral Method 

The Chebyshev spectral collocation method 15 is used to define the 
derivatives dfi/dx in equations (4). To define 3?/3x, we expand ? in terms of 
the Chebyshev polynomials 

T (x) = cos(n cos ^x) (28) 

n 

in the truncated Chebyshev series 

N 

F(x) = l F T (x) 
n==0 n n 

A 

where the F^ are the expansion coefficients of the series, 
x, we introduce the change of variables 

x = cos0, 0 < 9 < ir. (30) 

Putting equation (30) into (28) and introducing the resulting expression into 
(29) gives 

N A 

F(x) = 1 F cos n6, (31) 

n=0 n 


(29) 

To form a range on 


a Fourier cosine series. To discretize equation (31), we define a set of 
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collocation points Xj by 


7TJ 

X j = COS N > 


j 0,1,2, •••jN 


and the discrete form of (31) becomes 


F j = F(X j } = Jo F n COS N 


The inverse of (33) gives the F as 

n 


- 4 - I V 1 r cos SJJ 

NC n j=0 J 2 N 


where 


j = 0 or j = N 
1 < j < (N - 1) 


Examination of equations (33) and (34) shows that the can be efficiently 
evaluated using the fast Fourier transform. 20 

Next, we differentiate F in equation (33) with respect to x giving 


F'(x) = 


V F T'(x), 
L . n n 
n=l 


A form of equation (35) without derivatives of the Chebyshev polynomials is 
preferred, so we rewrite (35) in terms of another series 


F'(x) = 


J F (1) T (x) 
n n n 

n=0 


(36) 
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and then proceed to relate the coefficients of the two series. The following 
recursion relation exists between the Chebyshev polynomials and their 
derivatives 15 


T' T' 

n+1 _ n-1 = 2_ 

n+1 n-1 C n 

n 


(37) 


where 



n = 0 

. 

n > 1 


Putting equation (37) into (36) and algebraically manipulating the resulting 
expression gives 


F'(x) = 


N 


l 

n=l 



( 1 ) 

n-1 


T' - 
n 


; (i) 

l t 


N 


n=l 


2n 


n 


(38) 


Introducing equation (35) into (38) and simplifying then results in 


2n 


F = 
n 


, ; a) 

J n-1 n-1 


- F 


( 1 ) 

n+1 


(39) 


A (1) 

an expression for the F^ given the F^. The procedure for finding 
A ( 1 ) 

the F^ is initialized by setting 




<* / 1 \ A (l) 

and then solving for through Fq by back substitution. 15 Then, 

A (1) 

knowing all the F^ , the required spatial derivatives of F can be calculated 
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from equation (36). This procedure can again be done efficiently using the 
fast Fourier transform (FFT). 

In summary, the computational procedure to find 3?V3x is carried out as 
follows. First, given initial values of F. = F(x.), find the Chebyshev 
coefficients from equation (34) or the FFT. Next, compute the F v ' from 
equation (39). Then, knowing the F^ , compute 3F/3x from equation (36) or 
the FFT. Once 3F/3x and the source term H (see equations (4)) are known at 
t = 0, the solution may be advanced in time with an appropriate temporal 
integration scheme. The scheme developed for this work is discussed in the 
next section of the paper. 

Temporal Integrator 

Once values for 3F/3x and H are determined as described above, there 
remains a system of ordinary differential equations in time that must be 
solved for the dependent variable vector U. A number of algorithms were 
surveyed for integrating this system of 0DE"s including pure explicit schemes, 
pure implicit schemes, and mixed explicit-implicit schemes. The pure explicit 
schemes were in general unattractive because the stiffness of the ODE system 
made the algorithm inefficient. Pure implicit schemes were also precluded by 
the difficulty of developing spectral algorithms for the spatial derivatives 
evaluated implicitly. Hybrid explicit-implicit algorithms therefore appeared 
to offer the most attractive approach. Following References 5-7, (see page 2) 
the explicit-implicit split was formed by computing the source 
term H implicitly at the new time level, and computing 3F/3x explicitly at the 
old time level to allow application of the spectral approach. Having made 
this choice, the equations were then integrated in time using a two-stage 
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Runge-Kutta technique. The algorithm was developed as follows. 
We first discretize equation (4) as noted above giving 



= U i" At [®i + J + O(At) 2 

sp 


(40) 


where n is the old time level and n+1 is the new time level. The 
+n+l 

vector H is then expanded in a Taylor series in time. 


H n +1 


- H" + + 0(4t) 2 


or 


H n+1 = H n + K n (U n+1 - U n ) + 0( At) 2 


(41) 


where K n is the Jacobian of H, . Putting equation (41) into (4), 
simplifying the resulting equation, and then rewriting in delta form gives, 


[I + AtK n ]AU n+i = -At [ (4—)? + H n l 

L i J . l L '■9x ; i i J 

sp 


(42) 


where AIL = . Examination of equation (42) shows that the 

bracketed term on the left-hand side is a block-diagonal matrix, the blocks 
being n by n submatrices with n the number of equations in the system (4). 
Since the matrix in equation (42) is diagonal, equation (42) is the most 
easily solved for AU by inverting the blocks, i.e., 


AU^ +1 = -At [i + AtK^] _1 R^ 


(43) 
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where [ ] 


-1 


represents a block invert, 


and 


R 


n 


_ / 3F \n 


+ H 


sp 


(44) 


is the steady-state residual vector. To apply the two-stage Runge-Kutta 
technique, we divide equation (43) into a predictor and corrector step as 
follows: 


Predictor 


AU 


" +1 - -it [I + itK"]"- 1 R? 
1 L 1 J 1 


n-i-1 -.n 


U? +1 = u! 1 + AU? +1 
l i l 


(45) 


Corrector 

AU? +1 = -At [i + AtK?]" 1 R? +1 
1 L 1 J 1 

(46) 

U n+1 = U n + I f AU^ 1 + AU n+1 ). 
i i 2 v i l } 

Starting with initial conditions for J, equations (45) and (46) are used to 
advance the solution from time level n to time level n+1. The process is 
continued until steady-state conditions, defined as a reduction of ten orders 
of magnitude in the steady-state residuals, are reached* 

The magnitude of the time step in equations (45) and (46) is chosen based 
on the physical time scales present at any given time in the solution. The 
fluid-dynamic time step, At^, can be shown numerically to be limited by the 


Courant condition, 
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At 


Ax 

f j u | + c 


(47) 


The chemical relaxation time for a species i is given by^l 


t 

c 



Changes in this relaxation time are then given by 


(48) 


At 

c 


A(pf i ) 


(49) 


since w^ remains nearly constant over a time step* For accuracy, we require 
that the chemical time step be chosen such that no change in specific mass 
fraction greater than A occurs over that time step. Equation (49) then 
becomes 


At 

c 



(50) 


A is initially set at 0.0001 for the computations that follow. The 
computational time step At is then chosen to be the minimum overall grid 
points, of the fluid and chemical time step, i.e., 


At = min(At f , At^) (51) 

(We also note here, following our earlier discussion, that other authors 8 ^ 
restrict the truncation error or the growth of truncation error in the 
numerical solution to be below some set level, by the appropriate choice of 
time step. This is also an important constraint to achieve a given level of 
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accuracy. This constraint on At has not been applied in the present work, but 
such an approach should be included as the work further matures). 


Initial and Boundary Conditions 

The governing equations (4) are hyperbolic and require initial conditions 
at each point to start the calculation and boundary conditions at the inflow 
boundary. Initial conditions are computed by first specifying an inflow Mach 
number and estimating an outflow Mach number. The interior Mach number 
distribution is then assumed to have a spatial variation which is linear. The 
total pressure and total temperature are assumed to be constant throughout the 
domain. Finally, the initial flow is assumed to be isentropic, so that 
isentropic relations can be used to compute the static pressure and 
temperature; these conditions are found from 





(52) 


(53) 


Knowing the static temperature and pressure and Mach number, the velocity 
distribution can be computed, and the density distribution can be found from 
the equation of state. Since the inflow boundary flow remains supersonic, 
boundary conditions are specified thereby holding conditions fixed at their 


initial values 
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RESULTS 

Having now developed a numerical algorithm for solving the equations 
governing chemically reacting flows, that algorithm will now be used to model 
the reacting flow in a rapid expansion supersonic diffuser. A rapid expansion 
diffuser was chosen such that high concentration gradients existed near the 
inflow boundary, providing a formidable test for the method. Results from the 
present algorithm were compared with a benchmark calculation from an existing 
finite-difference chemical kinetics code to validate the method. The 
comparison also allowed a demonstration of performance of the high-order 
accurate spectral method on grids which were quite coarse compared to grids 
required in the finite-difference calculation. 

The rapid expansion diffuser is shown in figure 1. The diffuser is two 
units long, has an initial cross-sectional area of 0.79 and a final cross- 
sectional area of 3.14, The diffuser wall is defined, as noted, by a shifted 
sinusoidal. Flow is introduced to the diffuser at M = 1.4, a velocity of 1230 
m/s, a temperature of 1900 K, and a pressure of 0.081 MPa. The chemical 
composition of the inflow is defined to be a three-tenths stoichiometric 
mixture of hydrogen fuel and air. 

Starting from the initial state described above, the governing equations 
were solved, using the algorithm in a time consistent manner, until steady- 
state conditions were reached. Independent benchmark calculations were also 
carried out with a previously validated Adams-Moulton implicit finite- 
difference program. In the results which follow, comparisons between the two 
methods will be made, first to verify the new procedure, and second to 
demonstrate its high spatial resolution on relatively coarse grids. A 
comparison of methods showing a time history of the chemical species is given 
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in figures 2 through 4 for H 2 , O 2 , OH, and ^0, respectively. Results are 
presented at the first grid point interior to the inflow boundary where both 
the flow field and species gradients are a maximum. Agreement between the 
Runge-Kutta spectral and the finite-difference calculations is excellent in 
all cases. 

Next, we compare spatial results from the two methods once steady-state 
conditions have been reached. The finite-difference solution required 101 

grid points before a grid independent solution, defined as a graphically 
imperceptible difference in the steady-state result between the present grid 
and next coarser grid, was attained. Calculations using the Runge-Kutta 
spectral code were carried out on 17 and 9 point grids. Comparisons of 
steady-state results for the two methods are given in figures 5 through 10. 
Figure 5 shows a comparison of the axial velocity profiles in the diffuser. 
The 17 point spectral solution and the 101 point finite-difference solution 
agree quite well throughout the diffuser. The 9 point spectral solution 
slightly overpredicts the velocity near the inflow boundary, but agrees well 
throughout the remainder of the diffuser. The overprediction is likely due to 
the failure of the coarsest spectral grid to predict adequately the high 
gradients that exist at the beginning of the diffuser. Temperature 
comparisons, given in figure 6, follow similar trends, with the 17 point 
spectral solution agreeing well with the benchmark, and the 9 point solution 
also agreeing well, except near the inflow boundary. Identical trends also 
occur when axial pressure profiles are compared in figure 7. 

Comparisons of axial species distribution computed by the two methods are 
given in figures 8 through 10. Prediction of the H 2 mass fraction by the 
spectral method with 17 grid points agrees well with the finite-difference 
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solution throughout the diffuser as can be seen by examining figure 8. The 9 
point spectral solution underpredicts the H 2 mass fraction near the inflow 
boundary, again due to the high spatial gradient in fjj there, but agreement 
again becomes good away from the inflow boundary. The spatial distribution of 
0 2 mass fraction is given in figure 9. The gradients are not as large for 
this species since 0 2 is in excess, and both 17 and 9 point grids agree well 
with the finite-difference solution. The steady-state species distributions 
for OH and H 2 0 are given in figure 10. The spatial gradients are again high 
for both species near the inflow boundary, and trends similar to those for H 2 
are repeated here. Agreement is again quite good when comparing the 17 point 
spectral and finite-difference results. The 9 point spectral solution still 
underpredicts gradients near the inflow boundary, however. 

A final comparison of methods is given in figure 11 which shows the rate 
of reduction of steady-state residual with iteration count at the first 
interior grid point. Since the 17 point Runge— Kutta spectral and the 101 
point finite-difference calculations yield comparable accuracy and have the 
same minimum spatial step size, it is reasonable to assess the relative 
efficiency of the methods using the result given in this figure. Note that 
the residual reduction rate by the spectral code is significantly greater than 


that provided by the finite-difference code. The maximum residual (at any 
grid point) is reduced with the spectral code by ten orders of magnitude in 
only 2400 iterations, whereas the finite-difference code requires 6000 
iterations to achieve the same level of residual reduction. The more rapid 
rate of residual reduction translates directly into a superior overall 
computational efficiency for the spectral Runge-Kutta method. 



CONCLUDING REMARKS 


A numerical method has been developed for solving the equations governing 
chemically reacting flow fields. In this method, spatial derivatives are 
discretized using a Chebyshev spectral collocation technique. Species source 
terras are calculated using a global hydrogen-air finite-rate chemistry 
model. The resulting ordinary differential equations in time are integrated 
using a two-stage partial implicit Runge-Kutta scheme that is effective in 
handling the stiffness, due to chemical kinetics, that is present in these 
equations. A computer program has been written using the spectral Runge-Kutta 
algorithm, and this program has been used to compute the flow in a rapid 
expansion supersonic diffuser. The diffuser flow field was also computed 
using an existing implicit finite-difference reacting flow code that had been 
validated in prior analyses. Comparison of results from the two programs 
indicated that the Runge-Kutta spectral code was accurate in predicting both 
the time evolution of the dependent variables and the final steady-state 
results. In addition, the spectral algorithm produced equivalent accuracy on 
relatively coarse grids as compared to the fine grid required by the finite- 
difference calculations. Based on these results, it appears that the Runge- 
Kutta spectral method offers promise for improving our ability to model 
chemically reacting flow fields. 
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Figure 1. Rapid expansion supersonic diffuser test case 
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Figure 2. Comparison of time history of hydrogen mass fraction 




Figure 3. Comparison of time history of oxygen mass fraction. 
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Figure 4. Comparison of time histories of hydroxyl and water 
mass fractions. 
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Figure 5. Comparison of axial velocity profiles. 
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Figure 6. Comparison of axial temperature profiles. 




Figure 7. Comparison of axi 
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Figure 8. Comparison of axial hydrogen mass fraction profiles. 




Figure 9 


Comparison of axial oxygen mass fraction profiles. 
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Figure 10. Comparison of axial hydroxyl and water mass fraction 
profiles . 




Figure 11. Comparison of steady-state residual reduction rates of 
the methods . 
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A numerical algorithm has been developed for solving the equations describing 
chemically reacting supersonic flows. The algorithm employs a two-stage Runge-Kutta 
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finite-difference procedure with resultant significant gains in computational 
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